NofSteps=500000;

Transitions=[0.95,0.10;0.05,0.90];
disp(['The transition matrix is']);
disp(Transitions);

States=zeros(1,NofSteps+1);
ResidenceTimes_0=zeros(1,NofSteps+1);
ResidenceTimes_1=zeros(1,NofSteps+1);
ResidenceCount=zeros(1,2)';
ResidenceCount(1)=1;

for i=1:NofSteps
	switch(States(i))
		case {0}
			if(rand(1)<Transitions(1,1))
				States(i+1)=0;
				else
					States(i+1)=1;
			end
		case {1}
			if(rand(1)<Transitions(2,2))
				States(i+1)=1;
				else
					States(i+1)=0;
			end
	end
end

for i=1:NofSteps+1
	switch (States(i))
		case {0}
			if (ResidenceCount(1)==0)
				if(ResidenceCount(2)~=0)
					ResidenceTimes_1(i)=ResidenceCount(2);
				end
				ResidenceCount(2)=0;
				ResidenceCount(1)=ResidenceCount(1)+1;
				else
					ResidenceCount(1)=ResidenceCount(1)+1;
			end
		case {1}
			if (ResidenceCount(2)==0)
				if(ResidenceCount(1)~=0)
					ResidenceTimes_0(i)=ResidenceCount(1);
				end
				ResidenceCount(1)=0;
				ResidenceCount(2)=ResidenceCount(2)+1;
				else
					ResidenceCount(2)=ResidenceCount(2)+1;
			end
	end
end

ResidenceTimes_0(ResidenceTimes_0==0)=[];
ResidenceTimes_1(ResidenceTimes_1==0)=[];



HistIntv0=1:max(ResidenceTimes_0);
HistIntv1=1:max(ResidenceTimes_1);

figure;

subplot(2,1,1);
hist(ResidenceTimes_0,HistIntv0);
title('Histogram of residence time in state 0');
xlabel('Residence Time / times');
ylabel('Numbers');

subplot(2,1,2);
hist(ResidenceTimes_1,HistIntv1)
title('Histogram of residence time in state 1');
xlabel('Residence Time / times');
ylabel('Numbers');
